## Voting and racial attitudes (retired whites only)
## Brian T. Hamel and Bryan Wilcox-Archuleta
## First: 5 May 2019  
## Last: 19 March 2020

## Loading packages
## install.packages(c("estimatr", "tidyverse", "magrittr", "texreg", "gridExtra", "scales"))
library(estimatr)
library(tidyverse)
library(magrittr)
library(texreg)
library(gridExtra)
library(scales)

## Loading data
load("01_data/dta.RData")

## Create shell
shell = dta %>% filter(white == 1 & retired == 1) %$%
  expand.grid(racial_flux = seq(min(racial_flux, na.rm = TRUE), max(racial_flux, na.rm = TRUE), by = 1),
              pid7 = round(mean(pid7, na.rm = TRUE), digits = 2),
              ideo5 = round(mean(ideo5, na.rm = TRUE), digits = 2),
              female = round(mean(female, na.rm = TRUE), digits = 2),
              age = round(mean(age, na.rm = TRUE), digits = 2),
              faminc = round(mean(faminc, na.rm = TRUE), digit = 2), 
              educ = round(mean(educ, na.rm = TRUE), digits = 2), 
              pct_white = round(mean(pct_white, na.rm = TRUE), digits = 2), 
              pct_black = round(mean(pct_black, na.rm = TRUE), digits = 2),
              pct_unemployed = round(mean(pct_unemployed, na.rm = TRUE), digits = 2),
              pct_college = round(mean(pct_college, na.rm = TRUE), digits = 2),
              log_per_cap_inc = round(mean(log_per_cap_inc, na.rm = TRUE), digits = 2),
              gini = round(mean(gini, na.rm = TRUE), digits = 2),
              south = round(mean(south, na.rm = TRUE), digits = 2),
              non_rural = round(mean(non_rural, na.rm = TRUE), digits = 2),
              log_pop_density = round(mean(log_pop_density, na.rm = TRUE), digits = 2)) %>% 
  na.omit()

## Models and predicted probabilities
pres_dem = lm_robust(pres_dem ~ racial_flux + pid7 + ideo5 + female + age + faminc
                     + educ + pct_white + pct_black + pct_unemployed
                     + pct_college + log_per_cap_inc + gini + south + non_rural
                     + log_pop_density, data = dta %>% filter(white == 1 & retired == 1),
                     clusters = zipcode, se_type = "stata")

pred_pres_dem = cbind(predict(pres_dem, shell, 
                              se.fit = TRUE, type = "response"), 
                      shell)

house_dem = lm_robust(house_dem ~ racial_flux + pid7 + ideo5 + female + age + faminc
                      + educ + pct_white + pct_black + pct_unemployed
                      + pct_college + log_per_cap_inc + gini + south + non_rural
                      + log_pop_density, data = dta %>% filter(white == 1 & retired == 1),
                      clusters = zipcode, se_type = "stata")

pred_house_dem = cbind(predict(house_dem, shell, 
                               se.fit = TRUE, type = "response"), 
                       shell)

rr = lm_robust(mean_rr ~ racial_flux + pid7 + ideo5 + female + age + faminc
               + educ + pct_white + pct_black + pct_unemployed
               + pct_college + log_per_cap_inc + gini + south + non_rural
               + log_pop_density, data = dta %>% filter(white == 1 & retired == 1),
               clusters = zipcode, se_type = "stata")

pred_rr = cbind(predict(rr, shell, 
                        se.fit = TRUE, type = "response"), 
                shell)

affirm = lm_robust(affirm ~ racial_flux + pid7 + ideo5 + female + age + faminc
                   + educ + pct_white + pct_black + pct_unemployed
                   + pct_college + log_per_cap_inc + gini + south + non_rural
                   + log_pop_density, data = dta %>% filter(white == 1 & retired == 1),
                   clusters = zipcode, se_type = "stata")

pred_affirm = cbind(predict(affirm, shell, 
                            se.fit = TRUE, type = "response"), 
                    shell)

## Table of coefs., and save
##############
## TABLE A7 ##
##############
texreg(list(pres_dem, house_dem, rr, affirm),
       file = "03_output/retired.tex", 
       label = "retired",
       caption = "Racial Flux, Voting Behavior, and Racial Attitudes (Retired Whites)", 
       custom.model.names = c("\\textit{President}", "\\textit{U.S. House}", "\\textit{Racial Resentment}", 
                              "\\textit{Affirmative Action}"),
       custom.coef.names = c("Intercept", "Racial Flux", "Party ID",
                             "Ideology", "Female", "Age", "Family Income",
                             "Education", "% White", "% Black",
                             "% Unemployed", "% College", 
                             "log(Per Capita Income)", "Gini Coef.", "South",
                             "Non-Rural", "log(Pop. Density)"), 
       reorder.coef = c(2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1),
       custom.gof.names = c(NA, NA, "Observations", NA, NA),
       stars = c(0.05, 0.01, 0.001),
       digits = 3,
       center = TRUE,
       include.ci = FALSE,
       caption.above = TRUE,
       scalebox = 0.7)

## Plot, and save
pred_pres_dem = cbind(pred_pres_dem, outcome = "President")
pred_house_dem = cbind(pred_house_dem, outcome = "U.S. House")
pred_vote = bind_rows(pred_pres_dem, pred_house_dem) %>% 
  mutate(upper = fit + 1.96 * se.fit,
         lower = fit - 1.96 * se.fit)

vote_plot = ggplot(pred_vote, aes(x = racial_flux, y = fit, ymin = lower, ymax = upper)) +
  geom_line(color = "red4") +
  geom_ribbon(alpha = .2, fill = "red1") +
  facet_wrap(~ outcome, nrow = 1, scales = "free") +
  labs(y = "Pr(Vote Democrat)", 
       x = "") +
  geom_rug(data = dta %>% filter(white == 1 & retired == 1), 
           aes(x = racial_flux), inherit.aes = FALSE, sides = "b") +
  scale_y_continuous(labels = number_format(accuracy = 0.01)) +
  theme(legend.title = element_blank(),
        panel.spacing = unit(1, "lines"),
        axis.line.y = element_blank())

pred_rr = cbind(pred_rr, outcome = "Racial Resentment")
pred_affirm = cbind(pred_affirm, outcome = "Affirmative Action")
pred_att = bind_rows(pred_rr, pred_affirm) %>% 
  mutate(upper = fit + 1.96 * se.fit,
         lower = fit - 1.96 * se.fit)
pred_att$outcome = factor(pred_att$outcome, levels = c("Racial Resentment",
                                                       "Affirmative Action"))

###############
## FIGURE A3 ##
###############
att_plot = ggplot(pred_att, aes(x = racial_flux, y = fit, ymin = lower, ymax = upper)) +
  geom_line(color = "red4") +
  geom_ribbon(alpha = .2, fill = "red1") +
  facet_wrap(~ outcome, nrow = 1, scales = "free") +
  labs(y = "Predicted Attitude", 
       x = "Racial Flux") +
  geom_rug(data = dta %>% filter(white == 1 & retired == 1),
           aes(x = racial_flux), inherit.aes = FALSE, sides = "b") +
  scale_y_continuous(labels = number_format(accuracy = 0.01)) +
  theme(legend.title = element_blank(),
        panel.spacing = unit(1, "lines"),
        axis.line.y = element_blank())

main = grid.arrange(vote_plot, att_plot, ncol = 1, nrow = 2)
ggsave(main, file = "03_output/retired.png", height = 4, width = 4, units = "in", dpi = 600)

## Clear R 
rm(list = ls())